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Abstract 

We have studied the Jarzynski equahty (JE) in van der Pol and Rayleigh oscillators which 
are typical deterministic non-Hamiltonian models, but not expected to rigorously satisfy the JE 
because they are not microscopically reversible. Our simulations that calculate the contribution to 
the work W of an applied ramp force with a duration r show that the JE approximately holds for 
a fairly wide range of r including r — )• and r — )• oo, except for r ~ T where T denotes the period 
of relaxation oscillations in the limit cycle. The work distribution function (WDF) is shown to be 
non-Gaussian with the [/-shaped structure for a strong damping parameter. The r dependence of 
R (= - ksT In (e-l^^)) obtained by our simulations is semi-quantitatively elucidated with the use 
of a simple expression for limit-cycle oscillations, where the bracket (•) expresses an average over 
the WDF. The result obtained in self-excited oscillators is in contrast with the fact that the JE 
holds in the Nose-Hoover oscillator which also belongs to deterministic non-Hamiltonian models. 

PACS numbers: 05.70.-a, 05.45.-a, 05.40.-a 
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I. INTRODUCTION 



In the last decade, a significant progress has been made in theoretical studies on nonequi- 
librium statistics (for reviews, see Refs. [l-3|). The important three fluctuation theorems 



have been proposed: the Jarzynski equality (JE)^ 



the steady-state and transient fluctu- 



ation theorems [sHsl, and the Crooks theorem [7|, |8|. They may be applicable to nonequi- 
librium systems driven arbitrarily far from the equilibrium states. In this paper we pay our 
attention to the JE expressed by 
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(1) 



where W denotes a work made in a system when its parameter is changed, the bracket (•) 
expresses the average over the work distribution function (WDF), P{W), of a work performed 
by a prescribed protocol, AF stands for the free energy difference between the initial and 
final equilibrium states, and /3 (= l/ksT) is the inverse temperature of the initial state. 
Equation ([T]) includes the second law of thermodynamics, (W) > AF, where the equality 
holds only for the reversible process. The JE was originally proposed for classical isolated 
system and open system weakly coupled to baths which are described by the Hamiltonian 
[4] and the stochastic models |9|. Jarzynski later proved that the JE is valid for strongly 



coupled open systems 



10|. A validity of the JE has been confirmed by some experiments for 



systems which may be described by damped harmonic oscillator models H 



-l™. Stimulated 



by these experiments, many theoretical analyses 
with the use of the Markovian Langevin model 



lave been made for harmonic oscillators 



17 



20|, Fokker-Planck equation 2l|], and Hamiltonian model 



13 



16| , the non-Markovian Langevin model 



22 
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Recently the validity 



of the JE in nonlinear oscillators with anharmonic potentials has been investigated in Refs. 



18 



28|. 



In this paper, we will study the self-excited oscillators described by van der Pol and 
Rayleigh equations with state- and velocity-dependent dampings 29], 30 1. They are expressed 
by 



X 



-x-Cv + f{t), 
c(x^ — a) for the van der Pol model, 
clv"^ — a) for the Rayleigh model. 



(2) 
(3) 
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where f = x, a = 1, c (> 0) is a damping parameter, dot (■) denotes a derivative with 
respect to time t, ( stands for an auxihary variable, and f{t) expresses an apphed external 
force. Conditions of C > and C < express positive and negative dissipations, respectively. 
The van der Pol equation was proposed as a mathematical model of self-excited oscillations 
for a simple electric circuit with nonlinear triode valve The Rayleigh equation was 

introduced to show the appearance of sustained vibrations in acoustics 30||. Van der Pol and 
Rayleigh equations are formally equivalent in the sense that the van der Pol equation may be 
transformed to the Rayleigh equation and vice versa with a proper change of variables. These 
equations provide basic models for various nonlinear dynamics of systems in mechanical 
and electrical engineerings, biology, biochemistry and many other applications (for a recent 
review on nonlinear equations, see Ref. jsi])- Many studies have been reported on the van 



der 



^ol model, which may be regarded as a special case of the FitzHugh-Nagumo model 
331]. The properties of periodic solutions of the van der Pol oscillator are known in 



considerable detail for a sufficiently small or large damping coefficient. 

The van der Pol and Rayleigh oscillators belong to deterministic non-Hamiltonian models. 



The Nose-Hoover (NH) oscillator 3J, 1351] which has been widely adopted for a study of 



mo 

by 



ecu 
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ar dynamics also belongs to non-Hamiltonian models. The NH oscillator is described 



35] 
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-x-(v + f{t), 



(5) 
(6) 
(7) 



where C is a state variable of the thermal reservoir with the temperature T and tq stands 
for the relaxation time of (. It is noted that Eqs. ©-(HI) with a = ksT are similar to 
3qs. ©-(IT!) except for the fact that Eq. (jl]) is given by ^ while Eq. ([7]) is expressed by C, 
361 ]. The fluctuation theorem in non-Hamiltonian system coupled to thermostat has been 



discussed in Refs 



40|. Refs. 



41 



have provided the condition for the JE to hold 
in non-Hamiltonian (and Hamiltonian) model. The condition requires that the equilibrium 
canonical distribution should be given by j43| 



(8) 
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with 



ij{q) = -ksTlncPiq), (9) 

where H{Q) denotes energy of the system, 0(g) is the (normahzed) equihbrium distribution 
of bath variable q and A is an external parameter. The canonical distribution of the Nose- 



Hoover model with f{t) = X is given by 35 1 



P{X,V,0 OC e-'^(-'V2-A-W/2+rQCV2)^ (IQ) 



which satisfies the condition given by Eqs^ 



the JE holds in the Nose-Hoover model 



41 



8]) and (j9]) with Q = (x, v) and q = (, and then 
-|43|. In contrast, the equilibrium distribution of 



the van der Pol or Rayleigh oscillator which is an odd-shaped racetrack [44j, does not meet 
the condition given by Eqs. ([H]) and (|9]). Although this suggests that the JE does not hold 
in van der Pol and Rayleigh oscillators, it is worthwhile to examine how and to what extent 
the JE is violated in self-excited oscillators, which is the purpose of the present paper. 

The paper is organized as follows. In the next Sec. II, we briefly explain basic equations 
of van der Pol and Rayleigh oscillators. We examine the validity of the JE by simulations 
applying a ramp force with a duration r. Our simulations in Sec. Ill show that although the 
JE is not exactly satisfied in self-excited oscillators, the JE approximately holds in a fairly 
wide range of r values including r — (transient force) and t ^ oo (quasi-stationary force). 
Various types of analytical solutions for applied sinusoidal forces have been developed for van 
der Pol and Rayleigh models. It is, however, still difficult to obtain analytical solutions for 
arbitrary external forces including non-periodic ones. By using a simple analytic expression 
of solutions for an applied ramp force which is suggested by He's method for a limit cycle of 



self-excited oscillators j45|, we present in Sec. IV, a semi-quantitative analysis of the results 



of our simulations. Sec. V is devoted to our conclusion. 

II. SELF-EXCITED OSCILLATOR MODELS 
A. Energy, heat and work 

From Eqs. (l2])-(j4j), van der Pol and Rayleigh oscillators are described by 

x + x + Cv = f{t), (11) 



with 

{c(x^ — 1) for the van der Pol model, 
(12) 
c(f ^ — 1) for the Rayleigh model. 

When we set C = c (> 0) in Eq. ffTTj) . it expresses a damped harmonic oscillator. Multiplying 
X for the both sides of Eq. fill I) and integrating them over t, we obtain 

Uit) - f/(0) = Qit) + Wcit), (13) 

with 

Uit) ^ ^ + £f , (14) 

Q{t) = - [ Cx^dt, (15) 
Jo 

Wc{t) = [ f{t)xdt, (16) 



where U{t), Q(t) and Wc(t) stand for the internal energy, heat (dissipative energy) and 
classical work, respectively. Equation f lT3|) expresses the first law of thermodynamics. In 
order to show the JE, Jarzynski employed an alternative work defined by [4] 

Wj{t) = - [ jXt) x{t) dt, (17) 

which is related with Wc{t) as 

Wj{t) = -f{t)x{t) + /(O)x(O) + W,{t). (18) 

It is noted that U{t), Q(t), Wc(t) and Wj(t) depend on a microscopic history of the system 
of x{t) and v{t) for t > starting from their initial values of x(0) (= a;o) and f (0) (= Vq). 
Wj(t) has been employed for a study of the JE in this study. 

We have presented in the Appendix, some numerical calculations of thermodynamical 
quantities such as energy, heat and work of the van der Pol oscillator, which are evaluated 
both by single and multiple runs of simulations. It should be note that even for /(t) = 0, we 
obtain {U{t))o — {U{0))o 7^ in van der Pol (and Rayleigh) oscillators because of a dissipative 
contribution of {dQ(t)/dt)Q (see Figs. [T^ and in the Appendix), where {■)o stands for an 
average over initial states [Eqs. (1221) and (I23p ]. This is in contrast to the NH oscillator where 
the relations, (f/(t))o — (f/(0))o = and {dQ(t)/dt)o = 0, hold. This difference reflects on 
the difference in non-equilibrium properties of self-excited and NH oscillators: the JE does 
not hold in the former while it holds in the latter. 
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B. The Jarzynski equality 



For a study of the JE, we will apply a ramp force /(t) given by 



fit) = { 



for t < 0, 

g (^) for < t < r, 

g for t > r, 



(19) 



where r denotes a duration of the force and g its magnitude. By using Eqs. f lT7|) and f[T^ . 
we obtain a work induced by the applied ramp force, 



,y„ . WAr) = - (f) /' 



x{t) dt. 



(20) 



The WDF is expressed by 



P(W) = {S{W-Wo)), 



(21) 



where {■)o signifies the average over the canonically distributed {xq} and {fo}, 



satisfying the equi-partition relation given by 



(22) 



<^o>o = (^o>o = knT. 



(23) 



The JE in Eq. ([T]) may be rewritten as 



R = -iln(e-^^> = AF. 



/3 



(24) 



If the WDF is Gaussian, R is given by 



R = H 



2 ' 



(25) 



with 



a' = ((iy-/x)2). 



(26) 
(27) 



which stand for mean and variance, respectively, of the WDF. Of course Eq. (125|) is not 
valid for non-Gaussian WDF. 



III. MODEL CALCULATIONS 



A. The van der Pol oscillator 

Model calculations of the JE for van der Pol and Rayleigh oscillators will be reported in 
Sec. IIIA and IIIB, respectively. We have made simulations, solving Eq. f|TT]) by using the 
Runge-Kutta method with a time step of 0.0001 for initial states of {xq} and {vq} given by 
Eqs. ([22]) and ([23]) with ksT = 1.0. 

In order to get some insight into the van der Pol oscillator, we first show results without 
forces [f{t) = 0.0]. Time courses of x{t) and v{t) for a damping parameter of c = 10.0 
calculated by single runs with xq = 1.0 and vq = 0.0 are plotted in Figs. [T](a) and[T](b), 
respectively. Time courses of x{t) exhibit the relaxation oscillation with characteristic sharp 
periodic jumps. A parametric plot of x{t) vs. v{t) in Fig. [T](c) shows the limit cycle. The 
period of the relaxation oscillation depends on the magnitude of a damping parameter c. 
The dashed curve of Fig. [T](d) expresses the c dependence of period T for f{t) = 0.0, which 
is increased with increasing c. When a constant force / = 0.5 is applied to the oscillator, its 
oscillation period is further increased as shown by the solid curve in Fig. [11(d). 

Figure [2t^a) shows time courses of x{t) when a ramp force with r = 100.0 is applied to 
the van der Pol oscillator for c = 1.0 with initial conditions of Xq = 1.0 and Vq = 0.0. The 
period of the oscillation with the applied ramp force with g = 0.5 (solid curve) is gradually 
increased compared to that with g = 0.0 (dashed curve). Figure [2](b) shows a similar plot 
of x{t) for c = 10.0. The period of the oscillation with c = 10.0 is longer than that with 
c = 1.0 by a factor of about three. An applied ramp force with g = 0.5 and r = 100 induces 
a work Wj{t) whose time course is shown by the chain (solid) curve for c = 1.0 (c = 10.0) 
in Fig. Wic). We obtain Wq = -0.110 {Wq = -0.126) for c = 1.0 (c = 10.0) by a single run 
with initial values of Xq = 1.0 and Vq = 0.0. 

Calculating Wq in Eq. ( l20i) for given initial states, we have obtained the WDF, P{W), 
with the use of Eq. (I2T1) . whose results for r = 0.1, 1.0, 10.0 and 100.0 with c = 10.0 are 
plotted in Fig. [31 Although the WDF for r = 0.1 is Gaussian, those for r = 1.0, 10.0 
and 100.0 are non-Gaussian with the U-shaped structure, which are quite different from the 
Gaussian distributions obtained in harmonic oscillators. Figure [H shows P{W) for various 
values of c with a fixed r = 1.0. With decreasing the damping parameter from c = 10.0, 
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the WDF is changed from the double-peaked distribution to the single-peaked Gaussian-like 
distribution. 

Figures [S](a), [S](b) and[S](c) show r dependences of /i, a and -R, respectively, for c = 1.0 
(dashed curves), c = 5.0 (dotted curves) and c = 10.0 (solid curves) for ramp forces with 
g = 0.5. is almost zero for r = 0.1 and it gradually decreased to —0.125 for r = 1000.0. 
In contrast, a ^ 0.5 at r = 0.1 and it goes to zero at r = 1000.0 with a small bump at 
T ~ 3.0. The calculated R with c = 10.0 (solid curve) is in nearly agreement with AF 
(= —g^jl = —0.125) for r > 100.0 and r < 0.2, but it significantly deviates from AF for 
0-2 ^ T ^ 100.0. The discrepancy oi R ^ AF implies a violation of the JE. This deviation 
becomes less significant for smaller values of c = 5.0 and 1.0, and it vanishes for c = 0.0 
(harmonic oscillator) where the JE holds. 

B. The Rayleigh oscillator 

Next we study the case of the Rayleigh oscillator. Time courses of x{t) and v(t) of 
the relaxation oscillation for c = 10.0 and f(t) = 0.0 are plotted in Figs. [6](a) and[6](b), 
respectively, and a parametric plot of x{t) vs. v{t) in Fig. [6](c) exhibits the limit cycle. The 
period of the relaxation oscillation depends on the magnitude of c. The dashed curve of 
Fig. E](d) expresses the c dependence of the period T for f{t) = 0.0, which is increased with 
increasing c. The period of the oscillator for a constant / = 0.5 (solid curve) coincides with 
that for / = 0.0 (dashed curve) in Fig. [6](d). 

The dashed (solid) curve in Fig. El^a) shows x{t) of the Rayleigh oscillator with c = 1.0 
for ramp forces with g = 0.0 {g = 0.5) and r = 100.0, which is calculated by single runs 
with initial condition of Xq = 1.0 and Vq = 0.0. We note that x{t) for g = 0.5 is gradually 
shifted upward compared to that for g = 0.0. Figure Ulh) shows a similar plot of x{t) for 
the case of c = 10.0, whose period is longer than that for c = 1.0 in Fig. [Tl^a). The time 
course of Wj{t) for an applied ramp force with g = 0.5 and r = 100.0 is shown by the chain 
(solid) curve for c = 1.0 (c = 10.0) in Fig. [7](c). A work induced by the applied force is 
Wo = -0.123 {Wo = -0.101) for c = 1.0 (c = 10.0) for a given initial condition of xo = 1.0 
and Vo = 0.0. 

Calculated WDFs for various r are plotted in Fig. [HI WDFs for r = 10.0 and 100.0 have 
U-shaped structures, while they become the Gaussian-like distribution for r = 0.1 and 1.0. 
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Figures [9](a), [9](b) and[9](c) show r dependences of /x, o and i?, respectively, for c = 1.0 
(dashed curves), c = 5.0 (dotted curves) and c = 10.0 (sohd curves). Their r dependences 
are similar to those for the van der Pol oscillator shown in Figs. El^a), EJ^b) andO^c). We 
note that o for c = 10.0 has a large maximum at r ~ 10.0 where PiyV) has the two-peak 
structure as shown in Fig. [SI The calculated R for c = 10.0 (solid curve) is nearly in 
agreement with AF (= —0.125) for r > 20.0 and r < 0.5, but significantly deviates from 
AF for 0.5 ^ T < 20.0. This deviation is reduced for smaller c values of c = 5.0 (dotted 
curve) and 1.0 (dashed curve), and it vanishes for c = 0.0 which corresponds to a harmonic 
oscillator. 



IV. DISCUSSION 

A. WDF of harmonic oscillators 

We have tried to elucidate the result obtained by simulations having been reported in 
the preceding section. In a recent paper. He |45| has discussed a limit cycle of self-excited 



oscillators, by using a very simple expression for a relaxation oscillation given by 

x{t) = Acosut, (28) 

where an amplitude A and frequency u are determined as a function of a damping parameter 
c with the use of the variational method j^]- It has been shown that they are give i i by 
A = 2.0 and u = 3.8929/c + ©(c'^) for the van der Pol oscillator with c > 1 jsil, |45|. 



45| . we will study the WDF and the JE of self-excited oscillators 



Extending He's method 
in the following. 

Before discussing the WDF of self-excited oscillators, we briefly explain that of a harmonic 
oscillator [c = in Eq. f lTT]) ]. 



i + x = f{t), (29) 
whose solution for the applied ramp force f(t) is given by 



x{t) = Xo cos if: + f sin t + / sm{t — t')f{t') dt' , 

Jo 



(30) 



g(t — smt) , . 

Xo cost + fo smt H for < t < r. (31) 

r 



From Eq. (120!) . we obtain a work for a given initial condition of Xq and Vq, 



Wo = Cxo + Dvo + 



with 



C 



D 



g sm T 



g{l — cost) 
r 

, /(I -cost) 



2 r2 

The WDF in Eq. 1^ is given by 

P{W) = — / e*«w^(e-'"'^«)n 



(32) 

(33) 
(34) 
(35) 

(36) 



with 



exp 



iuCxt] 



dxo / exp 



— iuDvQ 



dvn 



oc exp(— ZM0jexp 



A simple manipulation with Eqs. (136|) and ( 137|) leads to the Gaussian WDF given by 



P{W) 



with 



V27ra2 

^2 ^2(^ _ pQg^^ 



2 2^2(^1 -cost) 



The average of e '^'^ over P(H^) in Eq. ([T]) is given by 



which yields 



^2 2 



(37) 

(38) 

(39) 
(40) 

(41) 

(42) 
(43) 



Equation (H3|) implies that the JE holds regardless of a value of r in harmonic oscillators 
13l-|27|. 



10 



B. WDF of self-excited oscillators 



We now calculate the WDF of self-excited oscillators. Taking into account Eqs. (128|) and 
flHT]) . we have assumed that the solution of the limit cycle in a self-excited oscillator is given 
by 

x{t) = A cos{iot -6) + ^(^-'''^t) for < t < r. (44) 

r 

Here A and u depend on c as in Ref. {45 1 and a phase 6 is determined by initial conditions 

of xq and vq, 

tane = ^, (45) 

UJXq 

which arises from 

Xo = Acos9, vq = uAsinO. (46) 

Note that A in a self-excited oscillator is assumed to depend on c but to be independent of 
initial condition of xq and fo, while A in a harmonic oscillator depends on them as given by 
A = ^xl + vl in Eq. (EI]). 

Substituting Eq. (jHj) into Eq. (l20l) . we obtain a work performed by the ramp force for 
a given ^, 

'^^\\- I Q\ ^ ■ Q\ ^2(l-cosr) 
— sm(wr -0) + sm 0\ - — ^ ^ . (47) 

The average of (e~*"'^°)o in Eq. fl36|) is given by 

/OO /"OO 
00 J —00 

Transforming Eq. ( I48l) to the polar coordinate and using Eq. (H6l) . we obtain 



with 



/i(^) = ( — ) [sin(wr -9)+ sin^], (50) 



WdCos{e~5), (51) 



cur 

(1 — coswr) , 

tan 5 = ^ (53) 

smwr 

, /(I -cost) 

= + • (54) 
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A substitution of Eq. (H9!) into Eq. (136!) leads to 



P{W) oc 



g-(MV2)(cos^ 



(55) 



which yields the WDF given by 



P{W) ~ 




1 



for -ly^ + /i < py < IVd + /i. 



(56) 



Equation (l56l) expresses the U-shaped WDF which is divergent at two edges of {±Wd + ii). 
Thus when c is increased from zero, the WDF changes from the Gaussian [Eq. fl55]) ] to 
U-shaped non-Gaussian [Eq. fl56l) ]. just as shown in Fig. HI 

Figures [TOlfa) and [TOT b) express and /x, respectively, as a function of r, which are 
calculated by Eqs. ([52]) and ([MD with A = 2.0, g = 0.5, u = 27t/T and T = 19.07 [Figs. 
Wl^d) and[6^d)]. We note in Fig. [TOT a) that the second law of thermodynamics holds because 
fi > AF (= —0.125). Wd has an interesting r dependence, which is similar to that of a in 
harmonic oscillator given by Eq. (I40l) . Two edges of {±Wd + /u) are plotted by solid curves 
in Fig. fTOT c). where squares (circles) express upper and lower edges of the WDF obtained 
by simulations for the van der Pol oscillator (Rayleigh oscillator). Oscillating behaviors in 
upper and lower edges of the WDF obtained in simulations are well reproduced in Fig. [TOT c). 

By using the WDF given by Eq. ( l56l) . we may evaluate (e~^'^) in Eq. ([T]), 



where In{z) expresses the modified Bessel function of the first kind. From Eq. ( 1571) . R in 
Eq. (I42l) is given by 



In the limit of r = oo, we have R = AF because fi = —g^jl = AF, Wd = 0.0 and 
/o(0) = 1.0. In the opposite limit of r = 0.0 where /z = 0.0 and Wd = gA, we obtain 
R = In Io{PgA) which is generally different from AF. By using Eqs. fl52|) . f l5^ and 

( 15S]) with A = 2.0, g = 0.5, co = 271 /T and T = 19.07, we have calculated R which is plotted 
by the solid curve in Fig. [TTJ For a comparison, we show by squares and circles, results 
of simulations for van der Pol and Rayleigh oscillators, respectively, with c = 10.0. We 
note that the r dependence of R for t > 1.0 obtained by simulations is semi-quantitatively 
explained by our analysis. 



{e-P"^) = e-^'^WWd) 



(57) 



R = fi--\nIo{f3Wd). 



(58) 
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Our calculation, however, yields poor results for r < 1.0 where x{t) in Eq. ( H^ is not 
a good approximation because an oscillation cannot become a limit cycle for t ~ r ^ T. 
Actually with decreasing r at r < 1.0, the WDF is changed from the U-shaped distribution 
to the Gaussian distribution, as shown in Figs. E] and [Hi Our calculation for r < 1.0 
may be improved if the WDF is phenomenologically interpolated between the Gaussian and 
U-shaped distributions as given by 

P^W) oc -jL= e-(^-'^)V2'^^ + ilzZ) ^ (59) 

with 

p = e-^l^\ (60) 

where tq denotes a parameter. The Gaussian and U-shaped WDFs are dominant for small 
and large r, respectively, and they are interpolated between small and large values of r with 
a factor p. The average of (e"^'^) is given by 

(e-/^^> = p e-^^^-'^^'/^) + (1 _ e-f'n^^ms). (61) 

The bold solid curve in Fig. [TT] expresses R obtained by Eqs. (H2|) . (160|) and (16T|) with 
To = 1.0. We note that i? deviates from AF for 1-0 < r < 10.0 although R ~ AF for 
r ^ 1.0 and r ^ 10, as shown by simulations. 



V. CONCLUSION 



Studying the JE in van der Pol and Rayleigh oscillators to which a ramp force with a 
duration r is applied, we have obtained the following results: 

(i) The JE nearly holds in a fairly wide range of r including transient (r — )■ 0) and quasi- 
stationary forces (r — )■ oo), although the JE is not rigorously satisfied |4l|-|43]. 

(ii) The WDF has the U-shaped structure for a large damping parameter, and 

(iii) The r dependence of R (= — /c^Tln <(e~'^^)) [Eq. fl24|) ] obtained by our simulations 
may be semi-quantitatively accounted for by our analysis with a simple expression of x(t) 
for a limit cycle whose amplitude is assumed to be determined by a damping parameter but 
not sensitive to initial conditions. 

The item fi) is in contrast with results of NH oscillators where JE holds 



41 



43|. Derivations 



of the JE require that the equilibrium canonical distribution in non-Hamiltonian systems 



13 



should satisfy the condition given by Eqs. ([8]) and (|9]) [43|. Van der Pol and Rayleigh 
oscillators do not meet the condition, while it is held in the NH oscillator 

lill-Q- Ahhough 

our simple analysis in the item (iii) may explain essential features of van der Pol and Rayleigh 
models, a development of more advanced theory is desirable for a better understanding of 
their properties. It would be interesting to examine our result by experiments, for example, 
by electrical circuits consisting of nonlinear elements. 
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Appendix: A. Energy, heat and work of the van der Pol oscillator 

In this Appendix we will present some model calculations of thermodynamical quantities 
such as the energy and heat in the van der Pol oscillator, which are evaluated both by single 
and multiple runs of simulations. Figures [T^ a). IT^ b). ITST c) and ITWd) show x(t), U{t), 
dQ{t)/dt and Wj{t), respectively, of the van der Pol oscillator with c = 1.0 for applied ramp 
forces with g = 0.0 (dashed curve) and g = 0.5 (solid curve) for r = 10.0 calculated by 
single runs with initial conditions of Xq = 1.0 and Vq = 1.0 [U{0) = 1.0]. The period of the 
oscillation with the applied ramp force with g = 0.5 is gradually increased compared to that 
with g = 0.0. We obtain U{t) - f/(0) = Q{t) for g = 0.0 where Wj{0) = Wc{t) = 0.0 in Eq. 
(IT^ . The heat (energy) flows from an environment to the oscillator for dQ{t)/dt > 0, and 
for dQ{t)/dt < the heat (energy) flow is reversed. Periodic energy exchanges are realized 
between the oscillator and environment in the limit cycle. For an applied ramp force with 
g = 0.5, Wj{t) is time dependent at < t < 10.0 and it becomes constant (= —0.191) at 
t > 10.0 where f{t) = 0, as shown in Fig. [TWd). 

Figures [T2]J^e)-[T2]J^h) show similar plots of relevant thermodynamical quantities in the van 
der Pol oscillator with a larger damping constant of c = 10.0 which are calculated also by 
single runs with the same initial condition of Xq = 1.0 and Vq = 1.0. The period of relaxation 
oscillation of x{t) for c = 10.0 is larger than that for c = 1.0. Time dependences of U{t) and 
dQ{t)/dt for c = 10.0 become much significant than those for c = 1.0: note that scales of 
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ordinates in Figs. and flW e) are much larger than those in Figs. [TW b) and [TW c). We 

obtain a positive Wj{t) (= 0.551) at t > 10.0 in Fig. MM)- 

Related thermodynamical quantities averaged over 100,000 runs with canonically dis- 
tributed initial states of {xq} and {vq} with ksT = 1.0 [Eqs. f l2^ and fl25]) ] are plotted 
in Fig. M Figures [ini^a), [ISl^b), CSl^c) and HSl^d) show (x(t))o, {U{t))o, {dQ{t)/dt)o and 
{Wj(t))o, respectively, of the van der Pol oscillator with c = 1.0 for applied forces oi g = 0.0 
(dashed curves) and g = 0.5 (solid curves) with r = 10.0. We note in Fig. [TST a) that 
although (x)o = for (7 = 0.0, (x)o for g = 0.5 at t > 10.0 expresses a small limit-cycle 
oscillation superposed on a constant of 0.5. This is because random initial states are effec- 
tively biased by an applied force. For g = 0.0, ([/(t))o = 1.0 at t = 0.0 and it becomes 
about 2.0 at t > 5.0 which is determined by amplitudes of x{t) and v{t) in the limit cycle, 
as shown in Fig. [TST b). It is noted that even for g = 0.0, we obtain (f/(t))o — (f/(0))o 7^ 0.0 
which is due to finite dissipative contributions of {dQ/dt)^ (7^ 0.0) between an oscillator and 
environment. Comparing Fig. [TST b) with Fig. [T2r b). we note that magnitudes of {U{t))o 
become much smaller than those of U{t) for a single run. Owing to an applied ramp force, 
{U{t))o for g = 0.5 is lower than that for g = 0.0. We obtain {Wj)o = -0.110 at t > 10.0 as 
shown in Fig. [TST d). 

Similar plots of relevant thermodynamical quantities for the van der Pol oscillator with a 
larger c = 10.0 are presented in Figs. [T3r e)-[T3rh). The initial averaged energy of (f/(t))o = 
1.0 at t = 0.0 is increased to about 2.4 ~ 3.6 at t > 5.0 for g = 0.0 in Fig. [THJ^f). This 
increase in {U(t))o is due to energy supplies from environment to the oscillator which are 
rapidly accomplished at < t < 1.0 both for g = 0.0 and 0.5 as shown in Figs. ITST f) and 
[T3](g). Figure [TST h) shows {Wj)o = —0.092 at t > 10.0, which is in contrast to a positive 
Wj = 0.551 for a single run shown in Fig. [T2r h). 



[1] C. Bustamante, J. Liphardt, and F. Ritort, Phys. Today 58, 43 (2005). 

[2] F. Ritort, Nonequilihrium fluctuations in small systems: From physics to biology^ in Advance 

in Chemical Physics, vol. 137, Ed. S. A. Rice (J. Wiley & Sons, Inc., 2008). 
[3] S. Ciliberto, S. Joubaud, A. Petrosvan. larXiv:1009.33621 
[4] C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997). 



15 



[5] D. J. Evans, E. G. D. Cohen, and G. P. Morriss, Phys. Rev. Lett. 71, 2401 (1993). 

[6] D. J. Evans and D. J. Searles, Phys. Rev. E 50, 1645 (1994). 

[7] O. Narayan and A. Dhar, J. Phys. A 37, 63 (2004). 

[8] G. E. Crooks, Phys. Rev. E 60, 2721 (1999). 

[9] C. Jarzynski, Phys. Rev. E 56, 5018 (1997). 

[10] C. Jarzynski, J. Stat. Mech.: Theor. Exp. P09005 (2004). 

[11] J. Liphardt, S. Dumont, S. Smith, I. Tinoco, C. Bustamante, Science 296, 1833 (2002). 

[12] G. M. Wang, J. C. Reid, D. M. Carberry, D. R. M. WiUiams, E. M. Sevick, and Denis J. 

Evans, Phys. Rev. E 71, 046142 (2005). 

[13] F. Douarche, S. Cihberto, A. Petrosyan, and I. Rabbiossi, Europhys. Lett. 70, 593 (2005). 

[14] F. Douarche, S. Joubaud, N. B. Garnier, A. Petrosyan, and S. Cihberto, Phys. Rev. Lett. 97, 

140603 (2006). 

[15] S. Joubaud, N. B. Garnier, F. Douarche, A. Petrosyan, and S. Cihberto, C. R. Physique 8, 
518 (2007). 

[16] S. Joubaud, N. B. Garnier, S. Cihberto, J. Stat. Mech., P09018 (2007). 

[17] F. Zamponi, F. Bonetto, L. F. Cughandolo and J. Kurchan, J. Stat. Mech. P09013 (2005). 

[18] T. Mai and A. Dhar, Phys. Rev. E 75, 061101 (2007). 

[19] T. Speck and U. Scifcrt, J. Stat. Mech. L09002 (2007). 

[20] T. Ohkuma and T. Ohta, J. Stat. Mech. PlOOlO (2007). 

[21] S. Chaudhury, D. Chatterjee and B. J Cherayil, J. Stat. Mech. P10006 (2008). 

[22] A. Dhar, Phys. Rev. E 71, 036126 (2005). 

[23] C. Jarzynski, Comptes Rendus Physique 8, 495 (2007). 

[24] C. Jarzynski, Eur. Phys. J. B. 64, 331 (2008). 

[25] R. Chakrabarti, arXiv: 0802.0268. 

[26] H. Hijar and J. M. O. de Zarate, Eur. J. Phys. 31, 1097 (2010). 

[27] H. Hasegawa, Phys. Rev. E 84, 011145 (2011). 

[28] A. Saha and J. K. Bhattacharjee, J. Phys. A 40, 13269 (2007). 

[29] B. van der Pol, Philos. Mag. Scr. 7 2, 978 (1926). 

[30] Lord Rayleigh, Philos. Mag. 15, 229 (1883). 

[31] Ji-Huan He, Intern. J. Mod. Phys. B 20, 1141 (2006). 

[32] R. FitzHugh, Biophys. J. 1, 445 (1961). 

16 



[33] J. Nagumo, S. Arimoto, and S. Yoshizawa, Proc. IRE 50, 2061 (1962). 
[34] S. Nose, Mol. Phys. 52, 255 (1984); J. Chem. Phys. 81, 511 (1984). 
[35] W. G. Hoover, Phys. Rev. A 31, 1695 (1985). 

[36] Equations ([2]), 1^ and ([Da) with a = UbT correspond to the equations proposed in Ref. [s^] 

for molecular dynamics with coupling to an external bath. 
[37] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak, J. 

Chem. Phys. 81, 3684 (1984). 
[38] G. Gallavotti and E. G. D. Cohen, Phys. Rev. Lett. 74, 2694 (1995). 
[39] D. J. Evans and D. J. Searles, Phys. Rev. E 52, 5839 (1995). 

[40] E.M. Sevick, R. Prabhakar, S. R. Williams, and D. J. Searles, Ann. Rev. Phys. Chem. 59, 
602 (2008). 

[41] M. A. Cuendet, Phys. Rev. Lett. 96, 120602 (2006). 
[42] M. A. Cuendet, J. Chem. Phys. 125, 144109 (2006). 
[43] E. Sch511-Paschinger and C. Dellago, J. Chem. Phys. 125, 054105 (2006). 
[44] B. L. HoUan, G. Ciccotti, W. G. Hoover, B. Moran, and H. A. Posch, Phys. Rev. A 39, 5414 
(1989). 

[45] Ji-Huan He, Phys. Rev. Lett. 90, 174301 (2003); 91, 199902(E) (2003). 



17 



FIG. 1: (Color online) (a) x{t), (b) v{t) and (c) a parametric plot of x{t) vs. v{t) in the van der 
Pol oscillator with c = 10.0 for f{t) = 0.0 with an initial condition of xq = 1.0 and vq = 0.0. (d) 
The c dependence of a period T with constant forces of f{t) = 0.0 (dashed curve) and f{t) = 0.5 
(solid curve). 



FIG. 2: (Color online) Time courses of x{t) of the van der Pol oscillator with (a) c = 1.0 and (b) 
c = 10.0 for ramp forces with g = 0.0 (dashed curve) and g = 0.5 (solid curve) for r = 100.0. (c) 
Wj{t) with c = 1.0 (chain curve) and c = 10.0 (solid curve) for g = 0.5 and r = 100.0. An applied 
ramp force f{t) is plotted by dotted curves in (a) and (b). Simulations are performed by single 
runs with initial conditions of xq = 1.0 and vq = 0.0. 



FIG. 4: (Color online) P{W) for c = 1.0, 2.0, 5.0 and 10.0 with r = 1.0 and g = 0.5 in the van 
der Pol oscillator. 



FIG. 5: (Color online) The r dependence of (a) (b) a and (c) R in the van der Pol oscillator with 
c = 1.0 (dashed curve), 5.0 (dotted curve) and 10.0 (solid curve) for ramp forces with g = 0.5, the 
arrow along the right ordinate in (c) expressing AF. The JE is expressed hy R = AF (= —0.125). 



FIG. 3: (Color online) P{W) for r = 0.1, 1.0, 10.0 and 100.0 with c = 10.0 and g = 0.5 in the 
van der Pol oscillator, P{W) for r = 10.0 and 100.0 being multiplied by factor of 1/2 and 1/20, 
respectively. 
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FIG. 6: (Color online) (a) x{t), (b) v{t) and (c) a parametric plot of vs. v{t) in the Rayleigh 
oscillator with c = 10.0 for f{t) = 0.0 with the initial condition of xq = 1.0 and vq = 0.0. (d) The 
c dependence of a period T with constant forces of f{t) = 0.0 (dashed curve) and f(t) = 0.5 (solid 
curve) . 

FIG. 7: (Color online) Time courses of x{t) of the Rayleigh oscillator with (a) c = 1.0 and (b) 
c = 10.0 for ramp forces with g = 0.0 (dashed curve) and g = 0.5 (solid curve) for r = 100.0. (c) 
Wj{t) with c = 1.0 (chain curve) and c = 10.0 (solid curve) for g = 0.5 with r = 100.0. An applied 
ramp force f{t) is plotted by dotted curves in (a) and (b). Simulations are performed by single 
runs with initial conditions of xq = 1-0 and = 0.0. 

FIG. 8: (Color online) P{W) for r = 0.1, 1.0, 10.0 and 100.0 with c = 10.0 and g = 0.5 in the 
Rayleigh oscillator, P{W) for r = 10.0 and 100.0 being multiplied by factors of 1/10 and 1/5, 
respectively, and that for r = 1.0 being shifted upward by 0.2. 

FIG. 9: (Color online) The r dependence of (a) fi, (b) a and (c) R in the Rayleigh oscillator with 
c = 1.0 (dashed curve), 5.0 (dotted curve) and 10.0 (solid curve) for ramp forces with g = 0.5, the 
arrow along the right ordinate in (c) expressing AF (= —0.125). 

FIG. 10: (Color online) (a) Wd and (b) ^ as a function of r calculated by Eqs. ()52p and (I54|) with 
A = 2.0, g = 0.5, u) = 2'k/T and T = 19.07. (c) The r dependence of (iPF^^ + ^) (solid curves) and 
that of upper and lower edges of the WDF obtained by simulations for van der Pol (squares) and 
Rayleigh oscillators (circles) with c = 10.0, dashed and chain curves being plotted only for a guide 
of eye (see text). 
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FIG. 11: (Color online) The r dependence of R calculated by Eq. (jSSh (solid curves) and Eq. 
()6ip (interpolation: bold solid curve) with A = 2.0, g = 0.5, uj = 2tt/T and T = 19.07, and those 
for van der Pol (squares) and Rayleigh oscillators (circles) with c = 10.0 obtained by simulations, 
dashed and chain curves being plotted only for a guide of eye. The arrow along the right ordinate 
expresses AF (= —0.125). 



FIG. 12: (Color online) (a) x{t), (b) U{t), (c) dQ{t)/dt and (d) Wj{t) in the van der Pol oscillator 
with c = 1.0, and (e) x{t), (f) U{t), (g) dQ{t)/dt and (h) Wj{t) with c = 10.0 for applied ramp 
forces with g = 0.0 (dashed curves) and g = 0.5 (solid curves) (r = 10.0) evaluated by single runs 
with the initial condition of xq = 1.0 and vq = 1.0 yielding f7(0) = 1.0. 



FIG. 13: (Color online) (a) (x(t))o, (b) (t/(t))o, (c) {dQ{t)/dt)o and (d) {Wj{t))o in the van der 
Pol oscillator with c = 1.0, and (e) {U{t))o, (f) {Q{t))o, (g) {dQ{t)/dt)o and (h) {Wj{t))o with 
c = 10.0 for applied ramp forces with g = 0.0 (dashed curves) and g = 0.5 (solid curves) (r = 10.0) 
averaged over 100,000 runs with ksT = 1.0 (= (C/(0))o). Results for g = 0.5 in (f) and (g) are 
shifted upward by five and ten, respectively, for a clarity of figures. 
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